IMPROVED  VELOCITY  MODELS  IN  WESTERN  CHINA  FOR 
TWO-  AND  THREE-DIMENSIONAL  FINITE  DIFFERENCE  MODELING 


L.  K.  Steck,  C.  R.  Bradley,  and  L.  E.  Jones 
Los  Alamos  National  Laboratory 

Sponsored  by  the  U.S.  Department  of  Energy 
Office  of  Nonproliferation  and  National  Security 
Office  of  Defense  Nuclear  Nonproliferation 
National  Nuclear  Security  Administration 

Contract  No.  W-7405-ENG-36 


ABSTRACT 


In  order  to  improve  our  understanding  of  the  tectonics  of  western  China,  we  are  developing  methods  for 
creating  improved  three-dimensional  (3-D)  lithospheric  structural  models  for  the  region.  Stratamodel®,  a 
computational  toolbox  for  modeling  fully  heterogeneous  rock  properties  in  3-D  and  Arclnfo,  a  widely 
accepted  surface  mapping  routine,  form  the  framework  within  which  we  build  our  models.  To  explore  the 
accuracy  of  our  new  velocity  structures,  we  compare  synthetic  seismograms  from  a  full-waveform  finite 
difference  algorithm  with  regional  data  available  from  the  IRIS-DMC  for  seismic  events  in  western  China. 

In  this  paper,  we  present  seismogram  modeling  results  from  an  earthquake  in  the  Xinjiang  province  of 
western  China,  demonstrate  improvements  to  previous  velocity  models,  and  describe  the  scheme  we  are 
using  to  create  2—  and  3-D  velocity  models  for  western  China. 

On  January  30,  1999,  at  03:51:07  GMT  an  earthquake  occurred  at  41.67  W,  88.46  N  near  the  Lop  Nor  test 
site  in  western  China.  Using  an  anelastic  finite  difference  code,  we  model  the  records  from  this  event  at  the 
stations  AAK,  MAKZ  and  WMQ.  Preliminary  two-dimensional  models  were  obtained  for  each  source- 
receiver  path  from  the  Cornell  Database,  and  these  are  fused  with  appropriate  one -dimensional  velocity 
models  from  published  receiver  function  work  or  our  own  surface  wave  inversions.  For  each  source- 
receiver  path,  the  Cornell  Database  provides  a  topography  profile,  a  depth-to-basement  (or  basin)  profile, 
and  a  Moho  profile.  To  date  we  have  not  included  topography  effects  in  our  modeling. 

The  Cornell  Database  contains  two  surfaces  that  define  the  interfaces  between  the  basin  and  crust,  and  the 
crust  and  Moho.  These  interfaces  form  the  primary  structural  divisions  in  the  Arc -Info  data  set.  The 
structural  model  is  then  improved  with  the  addition  of  geophysical  province  models.  The  specific  path  from 
source  to  receiver  is  then  forward  modeled  and  the  velocity  model  refined.  The  ultimate  goal  of  the  model 
maker  is  to  provide  a  tool  for  building  2-D  and  3-D  gridded  velocity  models  for  travel-time  and  full 
waveform  seismic  estimates. 

Given  that  any  particular  source  receiver  combination  may  cross  several  geophysical  provinces,  full 
waveform  modeling  is  incorporated  to  model  the  specific  path  effects  on  the  waveform.  We  have  modeled 
both  elastic  and  anelastic  realizations  of  the  crust  for  several  events,  and  present  a  memory -efficient 
scheme  for  modeling  the  effects  of  Q  in  three  dimensions. 
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OBJECTIVES 

One  of  the  fundamental  problems  facing  discrimination  and  location  of  earthquakes  and  explosions  is 
creating  accurate  structural  models  for  the  lithosphere.  Path  effects  on  the  seismic  waveform  can  obscure 
spectral  and  time  domain  discrimination  indices  and  poor  velocity  models  generate  ambiguity  in  event 
depth  and  location. 
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It  has  been  our  objective  over  the  last  year  to  generate  a  methodology  for  providing  improved  velocity 
models  for  western  China  and  to  test  this  methodology  via  the  method  of  finite  difference  waveform 
modeling.  Our  specific  objectives  are: 


1)  Provide  a  method  for  generating  improved  velocity  models  by  incorporation  of  1-D  and  2-D  velocity 
models  into  an  algorithm  to  merge  these  models  into  a  3-D  structure. 

2)  Generate  an  interface  for  extracting  these  models  into  a  useful  form  for  either  ray  tracing  or  full 
waveform  modeling. 

3)  Test  the  accuracy  of  the  models  via  finite  difference  and  travel  time  calculations  for  known  earthquake 
and  explosion  paths. 

4)  Perfect  a  memory -efficient  method  for  simulating  intrinsic  attenuation  (Q)  in  the  lithosphere  for  full 
waveform  simulation. 

To  date,  we  have  addressed  items  (1)  and  (2)  through  two  commercial  software  systems:  Stratamodel®,  a 
software  package  designed  to  incorporate  several  types  of  structural  information  for  the  creation  of  a  global 
3-D  lithospheric  model,  and  Arclnfo®,  a  data  base  algorithm  we  are  using  to  pass  lithologic  surfaces  to  the 
waveform  or  travel-time  modeler. 

We  have  selected  an  event  from  January  30,  1999,  that  occurred  near  Lop  Nor  in  western  China  to  test  our 
wave  propagation  algorithms  and  the  accuracy  of  the  velocity  model  maker.  The  results  of  our  analysis 
follow. 

RESEARCH  ACCOMPLISHED 


Finite  difference  modeling  of  the  1999  January  30  earthquake  near  Xinjiang  Province,  western 
China 

On  January  30,  1999,  at  03:51:07  GMT,  an  earthquake  occurred  at  41.67  W,  88.46  N.  Using  an  anelastic 
finite  difference  code,  we  model  the  records  from  this  event  at  stations  AAK,  MAKZ,  and  WMQ. 
Preliminary  2-D  models  were  obtained  for  each  source-receiver  profile  from  the  Cornell  Database,  and 
these  are  fused  with  appropriate  1-D  velocity  models  from  published  receiver  function  work  or  our  own 
surface  wave  inversions.  The  Cornell  Database  contains  two  3-D  structures  from  which  2-D  profiles  can  be 
obtained:  a  “Eurasia”  model  (EU)  and  the  University  of  Colorado  (CUB)  dataset.  The  latter  is  derived 
from  inversion  of  surface  wave  data  and  is  thus  appropriate  for  modeling  longer  period  energy.  Both  the 
EU  and  CUB  models  are  comprised  of  3-D  “topography”,  “depth  to  basement”  and  “depth  to  Moho” 
information.  To  date  we  have  not  included  topography  effects  in  our  modeling.  We  compare  results  from 
2-D  profiles  obtained  from  the  EU  and  CUB  to  synthetics  from  purely  one-dimensional  models.  For  each 
source-receiver  path  presented  here,  our  synthetics  comprise  a  bandwidth  from  roughly  0.02  Hz  to  1.0  Hz. 

Station  WMQ 

Station  WMQ,  at  roughly  247-krn  epicentral  distance,  is  located  on  the  northern  edge  of  the  Tarim  basin 
(Figure  1).  Seismograms  at  WMQ  were  modeled  using  the  CUB  and  the  EU  basement  and  Moho  profiles, 
along  with  several  simple  1-D  models.  The  primary  velocity  structure  used  for  this  source-receiver  profile 
is  based  on  a  7-layer  modified  Tien-Shan  model  (Kosarev  et  ah,  1993)  over  a  mantle  half-space.  Velocities 
for  a  1-D,  3-layer  average  Tien-Shan  model  were,  in  turn,  a  simplification  of  this  primary  velocity  model. 

Both  velocity  models  are  given  in  Table  I.  Several  alterations  to  the  7-layer  model  were  made,  including 
removing  the  crustal  low-velocity  zones  and  thinning  the  crust  by  about  10  km.  None  of  the  source/velocity 
model  combinations  tried  were  able  to  fit  both  the  body  and  surface  waves  simultaneously.  The  TW  and 
GER2a  models  fit  the  body  wave  arrival  times  and  polarities  well,  particularly  for  the  prominent  depth 
phase  and  Lg.  The  surface  waves  are  not  well  matched  with  this  source  model.  In  contrast,  the  CMT  and 
GER1  models  fit  the  surface  waves  well,  but  did  not  match  the  polarities  of  the  body  waves.  The  GER1, 
GER2,  and  CMT  depths  appear  to  be  too  deep  by  about  5  km,  based  on  depth  phase  arrival  time  and 
waveform  matching.  The  depth  phase  arrival  time  and  polarity  are  consistent  with  it  being  an  sP.  Figure  2 
compares  several  of  our  models  for  WMQ. 


Station  MAKZ 


We  model  station  MAKZ,  at  770-km  epicentral  distance,  with  a  similar  modified  Kosarev  structure  since 
its  source-receiver  path  passes  primarily  through  the  Tien  Shan  fold  and  thrust  belt  (Figure  1).  Crustal 
thickness  along  the  profile  varies  from  about  42  km  to  nearly  50  km  and  Moho  depth  increases  gently  from 
the  source  end  of  the  profile  to  the  receiver  end,  of  whether  the  complex  EU  or  simpler  CUB  basement  and 
Moho  are  used.  A  seven-layer  structure  based  on  the  Kosarev  model  (reference  above),  a  simplified 
average  three-layer  (basin,  basement,  and  mantle)  structure,  and  a  model  based  on  phase  velocity  inversion 
work  by  Curtis  and  Woodhouse  (1997)  were  tested. 

The  7-layer  Kosarev  model  (Table  I)  was  fused  with  depth  to  basement  and  Moho  information  from  the  EU 
and  CUB  datasets,  while  the  3-layer  ‘average’  model  was  tested  for  the  CUB  profile  only.  Results  for  the 
Kosarev  model  andCMT  solution  show  reasonably  good  fits  to  the  body  wave  data,  but,  as  we  observed  at 
station  WMQ,  surface  waves  appear  phase  shifted.  The  3-layer  average  Tien  Shan  model,  coupled  with  the 
CU  basement  and  Moho  structures,  produced  similar  results. 

The  modified  Curtis  and  Woodhouse  model  produced  surface  waves  that  were  obviously  fast  by  more  than 
10  seconds  for  both  the  TCW  and  GER  source  mechanisms.  Note  that  the  TCW  source  produces  body 
waves  of  incorrect  polarity  and  unusual  complexity  on  the  radial  component.  The  GER  source  produces 
correct  polarities  for  the  body  waves,  and  a  more  distinct  ‘beating’  on  both  the  vertical  and  radial 
component;  a  feature  observed  in  the  data. 

All  velocity  models  produce,  with  varying  degrees  of  success,  a  waveform  feature  in  the  Lg  phase  that 
resembles  'beating'.  The  presence  of  this  feature  in  the  synthetics  may  be  due  to  the  2-dimensional  nature  of 
the  modeled  Moho  structure. 


Table  I:  Velocity  Structures 


(Modified)  Kosarev  Model 

Kosarev  et  ah,  1993 

deDth 

Pvel 

Svel 

0.0 

4.19 

2.42 

2.44 

5.0 

5.83 

3.37 

2.72 

10.0 

5.60 

3.11 

2.72  LVZ 

20.0 

6.12 

3.54 

2.78 

30.0 

6.21 

3.59 

2.82 

40.0 

7.09 

4.10 

2.97 

Mantle 

7.80 

4.51 

3.37 

Average  Tienshan  Model 

0-5 

4.19 

2.42 

2.44 

5-50 

6.205 

3.587 

2.82 

Mantle 

7.61 

4.40 

3.28 

(Modified)  Curtis  and 

Woodhouse  Model 

Curtis  and  Woodhouse,  1998 

0-5 

5.93 

3.43 

2.768 

5-50 

6.17 

3.57 

2.88 

Mantle 

8.11 

4.69 

3.65 

Modeling  the  Lithosphere  of  Western  China 

Discrimination  and  location  of  earthquakes  and  explosions  depend  on  the  source  model  but  also  depend  on 
the  effects  to  the  seismic  waveform  as  it  propagates  through  the  crust.  The  problem  of  creating  accurate 
structural  models  for  the  crust  is  difficult  in  regions  where  data  are  sparse  or  non-existent.  We  have  chosen 
to  use  Stratamodel,  a  toolbox  created  by  Landmark  Graphics  Corporation,  to  model  the  3-D  lithospheric 
properties.  The  algorithm  incorporates  known  published  data  and  interpolates  or  extrapolates  properties 
based  on  a  set  of  arithmetic  or  geostatistical  rules. 

It  has  been  our  objective  over  the  last  year  to  generate  a  methodology  for  providing  improved  velocity 
models  for  western  China  and  to  test  this  methodology  via  the  method  of  finite  difference  waveform 
modeling.  Below  we  discuss  the  data  sets  we  are  using  and  the  method  for  their  incorporation  into 
Stratamodel  to  create  better  models  of  the  lithosphere. 

Stratamodel 

Stratamodel,  which  models  heterogeneous  rock  and  fluid  properties  in  three  dimensions  for  geological 
analysis  and  visualization,  was  originally  designed  to  model  oil  and  gas  reservoirs.  We  have  reformatted 
our  data  sets  for  western  China  to  mimic  Stratamodel’ s  required  format. 

There  are  three  types  of  models  required  by  Stratamodel  in  order  to  create  a  full  3-D  model  of  the 
lithosphere:  1)  the  Stratigraphic  Framework  Model,  2)  the  Well  Model  and  3)  the  Attribute  Model. 

Stratigraphic  Framework  Model 

Gridded  horizons,  layer  boundaries,  and  detailed  reservoir  characterizations  (typically  well  data)  form  the 
basis  for  Stratamodel’ s  3-D  modeling.  The  gridded  horizons  we  have  chosen  are  to  the  topography,  basin 
and  Moho  surfaces  accessible  via  the  Cornell  database  ('http://www.atlas.geo.cornell.edu/ima.html).  Figure 
3  shows  the  basin  and  Moho  relief  from  the  IPE  Cornell  data  set.  Previous  modeling  studies  of  path  effects 
through  western  China  (e.g.,  Jones  et.  al,  1998;  Bradley  and  Jones,  1998)  have  used  these  surfaces  as  a 
framework  around  which  individual  crustal  models  (e.g.  Romanowitz,  1982;  Kosarev  et  al.,  1993;  Curtis 
and  Woodhouse,  1997;  Jih,  1998;  Mahdi  and  Pavlis,  1998)  have  been  inserted  in  a  “layercake”  manner  for 
the  particular  region  of  interest. 

These  data  surfaces  divide  the  lithosphere  into  sequences  that  have  undergone  distinct  geophysical 
processes  in  time,  pressure,  and  temperature  and  logically  should  be  treated  differently  by  the  model 
creation  process.  These  surfaces  are  the  reference  for  the  stratigraphic  framework  in  the  lithospheric 
model.  These  surfaces  dictate  the  cell  resolution  of  the  initial  model  and  subdivide  the  model  into 
geophysically  independent  sequences.  These  sequences  are  the  basin  (and  topography)  structure,  crustal 
structure,  and  Moho  or  lower  lithosphere  structure.  The  “well”  model  controls  layering  within  each  of 
these  sequences. 

Well  Model 

In  an  oil  reservoir  there  may  typically  be  well  data  to  enhance  the  information  on  the  geostratigraphy.  In  a 
similar  vein,  there  are  many  detailed  studies  of  regions  of  interest  in  western  China,  some  encompassing 
extensive  regions  (Jih,  1998;  Li  and  Mooney,  1998;  Mooney,  1999).  These  individual  studies  provide  1-D 
and  2-D  information  about  specific  geophysical  provinces.  Figure  4  shows  the  geophysical  provinces  we 
have  chosen  for  the  initial  model  of  western  China.  There  are  33  distinct  provinces,  each  of  which  has  a  1- 
D  structural  model  associated  with  it.  The  preliminary  velocity  models  were  derived  from  Jih,  1998,  and  Li 
and  Mooney,  1998.  Currently  there  are  fifteen  1-D  models  for  the  33  provinces  so  there  is  some 
redundancy,  but  we  have  allowed  for  expansion  when  our  knowledge  of  the  crust  improves. 

The  1-D  structural  models  are  treated  in  Stratamodel  like  wells.  They  provide  the  detailed  information  for 
the  sub-layering  within  each  sequence.  This  information  is  called  the  “attribute”  model. 


Attribute  Model 

Once  the  stratigraphic  framework  is  built  to  describe  structure  and  stratigraphy  and  well  models  is  in  place 
for  the  correlation  of  individual  provinces,  the  attribute  model  must  be  defined.  The  attribute  model 
describes  how  the  geophysical  rock  parameters  are  distributed  throughout  the  model  grid.  This  distribution 
can  be  based  on  simple  mathematical  relations  or  higher  order  geostatistical  parameters. 

Within  each  model  cell,  the  search  radius,  the  mathematical  bias  relative  to  the  “well”  data  and  the 
interpolation  scheme,  determines  the  attributes  associated  with  it.  For  our  current  model,  we  are 
calculating  Vp,  Vs,  density,  and  the  intrinsic  attenuation  coefficients  for  P-waves,  Qp,  and  S-waves,  Qs. 

The  Model 

The  stratigraphic  framework  defined  by  the  Cornell  surfaces  and  the  well  data  defined  within  each 
geophysical  province  have  given  an  initial  lithospheric  model  for  China.  Each  sequence  defined  by  the 
topography,  basin  and  Moho  surfaces  has  been  subdivided  by  four  layers  (Figure  5).  The  number  of  layers 
in  the  crustal  sequence  was  chosen  to  match  the  models  presented  in  Mooney,  1999.  The  basin  and  lower 
lithosphere  layering  will  be  improved  as  more  data  are  incorporated  into  the  model.  In  Figure  5  the 
topography  and  two  orthogonal  cross  sections  through  the  model  are  displayed.  In  the  basins,  the  layering 
is  assumed  to  be  proportional.  This  forces  the  layers  within  the  basin  to  pinch  out  as  the  basins  shallow. 
Within  the  crust  and  lower  crust,  the  layers  follow  the  top  of  the  sequences  that  define  them.  The  resulting 
model  follows  structural  trends  more  closely  than  simple  “layercake”  models. 

The  well  data  are  shown  in  Figure  6  for  the  33  provinces  used  in  this  model.  The  effect  of  the  wells  is  to 
bias  the  nearby  model  cells  with  information  from  the  physical  province.  Each  sequence  is  independent  of 
the  others.  This  forces  the  cells  that  lie  within  a  sequence  to  use  only  well  data  lying  between  the  surfaces 
that  define  that  sequence.  Currently,  the  weighting  scheme  for  the  cell  attribute  is  a  simple  inverse  distance 
weighted  average  of  the  well  data  and  a  search  radius  from  each  cell  equal  to  half  the  total  model  size. 

Future  models  will  refine  the  weighting  scheme  and  well  density  to  more  strictly  enforce  the  transition 
from  province  to  province  in  the  geophysical  model.  With  the  current  model  we  are  performing  full 
waveform  modeling  and  travel-time  tests  to  determine  the  accuracy  of  the  model. 

CONCLUSIONS  AND  RECOMMENDATIONS 


We  have  performed  2-D  modeling  of  an  event  near  Xinjiang,  China  at  stations  AAK,  MAK  and  WMQ. 
Depending  on  our  choice  of  focal  mechanism,  we  have  achieved  a  good  fit  for  the  body  wave  phases,  event 
depth  and  surface  waves.  No  single  mechanism  fits  all  of  these  features  simultaneously.  We  conclude  that 
path  effects  on  the  recorded  waveform  are  not  yet  being  fully  simulated  and  that  better  crustal  models  need 
to  be  developed.  These  are  our  recommendations. 

1)  A  robust  model-making  method  is  critical  to  address  the  need  for  2-D  and  3-D  structure  in  the  simulation 
of  path  effects  on  the  recorded  waveform. 

2)  Validation  through  modeling  of  the  velocity  structure  should  include  full  3-D  simulation,  and  there  may 
be  a  need  for  intrinsic  attenuation  in  the  simulations. 
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Figure  1 :  Map  of  source  region  for  the  30  January  1999  event.  Focal  mechanisms  are  shown 
as  predominantly  thrust  component.  Stations  modeled  are  indicated  by  the  back  azimuth  lines 
shown. 
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Figure  2. This  shows  (a)  vertical,  (b)  radial,  and  (c)  long-period  vertical  synthetic 
waveforms  for  several  source/velocity  model  combinations.  Instrument-corrected 
velocity  data  are  shown  with  the  synthetics. The  source  model  and  velocity  model 
are  indicated,  along  with  the  type  (if  any)  of  depth-to-basement  and  depth-to-moho 
information. The  'Thin'  velocity  model  is  an  abbreviated  flat-layer  Kosarev  model  with 
a  40-km-thick  crust  (no  E  U  or  C  U  B  information  used).  The  '3-layer'  model  is  an  average 
Kosarev  model  with  3  flat  layers. 


Figure  3.  The  IPE  basement  and  Moho  surfaces  from  Cornell.  These  surfaces  are  used 
to  define  the  initial  lithospheric  sequences.  Layering  within  the  sequences  is  dictated  by 
1-D  models  describing  distinct  geophysical  provinces. 
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Figure  4.  Defense  Advanced  Research  Project  Agency  (DARPA)  boundaries  used 
in  the  initial  velocity  model.  These  33  geophysical  provinces  have  a  1-D  velocity  model 
associated  (not  necessarily  unique)  and  are  used  to  enforce  a  structural  change  in  the 
model  generated  by  Stratamodel. 


Figure  5.  Velocity  model  for  China  from  Latitude  22  N  to  56  N,  Longitude  65  E  to  130  E . 
Topography  shown  as  the  uppermost  surface.  Two  orthogonal  cross  sections  show  the 
depth  and  detail  of  the  current  model:  12  heterogeneous  layes  in  3  sequences,  depth 
extent  125  km. 


Figure  6.  The  topography  has  been  stripped  away  from  figure  4  to  show  the  layer  sequences 
in  more  detail.  The  33  1-D  velocity  models  are  shown  in  as  "wells". 


